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Abstract 

Spectral methods are an efficient way to solve partial differential equations 
on domains possessing certain symmetries. The utility of a method depends 
strongly on the choice of spectral basis. In this paper we describe a set of bases 
built out of Jacobi polynomials, and associated operators for solving scalar, 
vector, and tensor partial differential equations in polar coordinates on a unit 
disk. By construction, the bases satisfy regularity conditions at r = 0 for 
any tensorial field. The coordinate singularity in a disk is a prototypical case 
for many coordinate singularities. The work presented here extends to other 
geometries. The operators represent covariant derivatives, multiplication by az- 
imuthally symmetric functions, and the tensorial relationship between fields. 
These arise naturally from relations between classical orthogonal polynomials, 
and form a Heisenberg algebra. Other past work uses more specific polynomial 
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bases for solving equations in polar coordinates. The main innovation in this 
paper is to use a larger set of possible bases to achieve maximum bandedness of 
linear operations. We provide a series of applications of the methods, illustrating 
their ease-of-use and accuracy. 
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1. Introduction 

Cylindrical polar coordinates find applications in countless areas of science and 
engineering. Important applications include pipe flow, laboratory studies of 
thermal convection, astrophysical accretion disks, electromagnetic waveguides, 
elastic deformation of rods, astronomical instrumentation, and plasma toka- 
maks. Many applications require the accurate and efficient solution of systems 
of partial differential equations (PDEs). Pseudospectral methods of different 
types prove useful for this task in many different geometries. In polar coordi¬ 
nates, the periodic nature of the azimuth angle allows the effective use of Fourier 
series, where 



( 1 ) 


After the Fourier transform, differentiation in 9 becomes multiplication by im. 


While Fourier analysis easily dispatches the azimuthal coordinate for functions 


on a disk, the radial coordinate presents difficulty for the following reason. For 
functions analytic everywhere on the disk, including the origin. 


fm(r) ~ as r 0, 


( 2 ) 


where F(r‘^) is an even function of r that is analytic at the origin. 

The coordinate singularity at the disk centre requires an m-th order zero for in¬ 
finite differentiability [SSlIll. Enforcing this condition in numerical calculations 
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presents challenges; especially for m ^ 1. Many authors address this chal¬ 
lenge with equally as many different techniques. Even considering regularity 
at the origin, the disk geometry allows a large number of possible orthogonal- 
polynomial bases mm- Zernike (1934) [32] produced the first practical set of 
polynomials for expanding functions on the unit disk. This basis proves partic¬ 
ularly useful in optical applications. Bhatia & Wolf (1954) [T| pointed out that 
this set is the only out of a possible infinity that contains “simple properties 
strictly analogous to that of Legendre polynomials. ” 

Boyd & Yu (2011) |3] provide a comprehensive review of the history and con¬ 
temporary methods used to solve Poisson’s equation in a disk. In particular, the 
paper reviews bases using Zernike-type polynomials, as well as the more common 
Chebyshev polynomials. The results for Chebyshev series range from accept¬ 
able to untenable. The diversity of Chebyshev methods results from different 
ways to represent the pole condition and/or the reflectional symmetry near the 
origin. A minimalist approach happens to produce the best option. This option 
expands even/odd-m modes in terms of an even/odd-degree Chebyshev series. 
This approach double wraps the disk using a Chebyshev series over — 1 < r < 1 
[31 [TO]. Compared to other Chebyshev options, simple even-odd matching 
works well with no other special intervention m- Even-odd matching and/or 
double-covering can satisfy equation ([^ with good-to-moderate accuracy. These 
schemes however do not enforce the analytic condition explicitly. This implies 
that singularities can still arise in higher-order derivatives; also see m- Even 
weak singularities can produce instabilities at the origin when performing time- 
evolution simulations. As a third option, the Roberts basis combines an even 
Chebyshev series with an explicit r™ prefactor. In spite of initial attractiveness 
(e.g., possessing a fast transform), this basis suffers from extreme numerical ill 
conditioning, and is not recommended [ 3 ]. 

Regarding the Zernike-type bases, Boyd & Yu point out that they are “More 
accurate for large m.” They also discuss the less-fortunate fact that Zernike 
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bases do not admit a fast transform in the radial direction. But that for various 
reasons “the advantages of ‘FFT-ability’ is not huge.”. They conclude that “It 
is difficult to definitely endorse one particular method for the disk because of 
the vast diversity of solutions to interesting engineering and science problems.” 
Furthermore, Slevinsky (2016) has recently made significant progress toward 
designing an effective fast transform from values at Chebyshev points to Jacobi 
coefficients that would work with the Zernike basis [3D]. For these reasons and 
more, we believe that polynomial bases that satisfy equation ([^ are very useful 
in many applications and are worthy of more detailed understanding. 

In addition to scalar-valued functions, many situations also require vector and 
tensor fields. Vectors introduce additional complications near the coordinate 
singularity. Much less work exists addressing these issues. In particular, the 
TO-th Fourier components of a vector field behave such that 

r ^ 0, (3) 

where, like F(r^) in equation ([^, V(r^) is an even function of r that is ana¬ 
lytic at the origin. We can (for example) see the necessity of equation (H by 
differentiating equation ([^ with respect to r. 

Sakai & Redekopp (2009) |3T] circumvented this issue by working with rescaled 
variables of the form which behaves like equation ([^. Li, Livermore 

& Jackson (2010) [16] use a poloidal-toroidal formulation to create a genuine 
(higher-order) scalar system out of a specific vector system. Using a technique 
equivalent to the r rescaling, Matsushima & Marcus (1995) [IS] show (and Boyd 
& Yu [3] reiterate) that Zernike polynomials produce pentadiagonal matrices 
for the solution of the radial portion of Poisson’s equation. Lastly, Townsend, 
Wilber & Wright (2016) develop an efficient low-rank approximation of scalar 
and vector functions on the disk that preserve regularity [33]. These methods 
work well for data analysis. Their application to time evolving systems remains 
less clear. 


4 


We show in this paper the non-necessity of radial rescaling and/or equation re¬ 
formulation. Previous works found recursion relationships for elementary oper¬ 
ators and rd/dr. Our calculus finds simpler factorisations of these operations 
in terms of r, d/dr and m/r. These are all of the elementary operators needed 
for full tensor calculus. This not only makes calculations easier to formulate, 
but also more numerically efficient and stable. In the process, we also show 
how to construct solutions to Poisson’s equation on the disk using only tiidi- 
agonal (as opposed to pentadiagonal) matrices; see the discussion in example 1 
in for more details. The foundation of the new results rests on exploiting a 
more general class of orthogonal polynomials. That is, we choose different bases 
to represent domain and range spaces of operators, so that coupling becomes 
banded. This mirrors using ultra-spherical polynomials for solving equations 
on the unit interval [3 [a [23]. Moreover, we incorporate azimuthally symmetric 
variable coefficients without destroying bandedness. This occurs via approxi¬ 
mating non-constant coefficients with finite-degree polynomials similar to |23j . 

A central theme of this paper demonstrates that increasing the collection of 
available bases can increase (i) the simplicity of a calculation’s numerical im¬ 
plementation; (ii) the speed to compute a solution; and (Hi) the accuracy of the 
result. We outline our following results: §2 derives properties for useful bases 
for polar coordinates (using properties of Jacobi polynomials). §3 shows how 
these bases respond to the covariant derivative operator in polar coordinates. 
§4 discusses multiplication by radial functions. §5 shows how the different bases 
relate to each other, and how the different operators form a Heisenberg Lie alge¬ 
bra. §2-5 build the necessary tools to represent and manipulate scalars, vectors, 
and tensors in different calculations. §6 applies these tools to a series of four 
example problems. §7 gives concluding remarks. 
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2. Calculus of Jacobi polynomials 


Throughout this paper, we use the definition 


z = 2r^- 1, 


( 4 ) 


where 0 < r < 1, implies —l<z<+l. The two following formulae provide 
the foundation for a simple calculus of functions on the unit disk. For any 
differentiable function F{z), 


d m 
dr r 

d 

dr r 


r^F{z) = ir^+^—Fiz), 

dz 


r^Fiz) = 2r 


m— 1 


TO + (1 + z) 


dz 


F{z). 


(5) 

( 6 ) 


The z-differential operators on the right-hand side of equations act on 

Jacobi polynomials in a particularly simple way. 


Jacobi polynomials P^’^\z) are the two-parameter set of classical orthogonal 
polynomials on —1 < z < 1 under the weight (1 — z)“(l -I- z)*”, where —1 < a, b: 

y''^P^^)(z)p(^^)(z)(l-zr(l + z)'’dz = (7) 

and 

^.6 ^ 2 °+'’+! r(n + a + l)r(n + 6+l) 

" 2n-|-a-|-6-|-l r(n -I- a -I- 6 -I-1) n! ’ 

where Jn,n' represents the Kronecker delta. Jacobi polynomials satisfy a Sturm- 
Liouville differential equation, and a three-term recursion relation. They also 
satisfy many other useful formula^ [55] . For the purposes of this paper we use: 


For differential operators, 


b+{l + z) 





(n + a + &+l)Pt+'’'’+'Hz), (9) 

(n + 6)p(“+i-''-i)(z). (10) 


^All relevant formulae are found in Digital Library of Mathematical Functions, Chapter 
on Orthogonal Polynomials, §18.3—§18.9: http://dlinf.nist.gov/18 
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For algebraic operators, 


(l + ^)P^^)(z) 


2(n+l) p(o.fc—1)/ 'I I 2(rt+b) p(a,fc— 

2n+a+h+l n+1 "C 2„+a_|_t,+ i ^ n 


( 11 ) 




n+g+b+l p{a+l,b)( 

2n+a+b+l n . 


n+b p(“+l.&) ( \ 
2n+a+b+l-^n-l 


( 12 ) 


For restriction operators, 


pw>(* = i)=(“+”)^ 


(13) 


Exchanging the roles of a and b, we can produce similar relations as equa¬ 
tions © -([l^, using 




(14) 


Doha & Bhrawy (2006) advocate for using Jacobi polynomials of varying degree 
to solve a number of problems in one and two dimensions. They specifically state 
that “The key to the efhciency of these algorithms is to construct appropriate 
base functions, which lead to systems with specially structured matrices that 
can be efEciently inverted.”^. We expand on this philosophy and construct a 
hierarchy of basis functions that conform particularly well to vector calculus in 
circular geometry. 


We define the following radial basis elements in terms of Jacobi polynomials. 




p' 
I rj 




(^). 


\/i^ 






22+fc+m ’ 


which implies an orthonormal basis 


^0 


(15) 


(16) 


From now on, the m index corresponds to the Fourier mode number in equa¬ 
tion Q. Factorials replace Gamma functions in equation (151 for integer values 
of k. Each element of the radial basis satisfies the analytic property in equa¬ 
tion ([^. We could consider a non-integer shift in the starting value of m in 
equation e.g., this could correspond to starting with Worland polynomi¬ 

als rather than Zernike m- The following analysis would work with such an 
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alteration, albeit with a large sacrifice to simplicity. We therefore only discuss 
integer-m values. 


The derivative operators for Jacobi polynomials imply the derivative operators 
on the radial basis 


d m 
dr r 
d m 
dr r 


^ 2y/n{n + k + m + l) , (17) 

Qt’^ir) = 2V(n + m)(u + fc + l)Q^i’™-i(r). (18) 


Equations 0 & ( [T^ form the foundation of a compact calculus of fields 
panded in terms of (5^’™(r). 


ex- 


Translating explicit statements like equations ([T7| & (181 into more abstract 


operators acting on vector spaces allows easier interpretation of many of the 
results and derivations. Therefore, we use Dirac’s bra-ket notatioij^ for the 
following analysis [3]. The following row-vector definition forms the foundation 
for our algebra, 


(fc. 


m.r \ = 


QS’™(r),gt’™(r), ... 


(19) 


The k, m indices parameterise different bases. The corresponding ket is the 


transpose | A:, m, r ^ ^ fc, to, r | .In the bra-ket notation equations (17) & (181 

become 


1 

71 


d TO 
dr r 


(fc, TO,r| = ( fc-I-1, TO ± 1, r 


( 20 ) 


where represent single-band matrices. The operators depend on the 
k,m indices, but from now on we omit these labels when the context makes 
their values clear. D~ is a diagonal matrix, and contains only a single band 
on the first super diagonal. The special case to = 0 permits only the raising 


operation, and D = . In equation (20), we choose the convention that 


r-space actions operate on the left, and matrix actions operate from the right. 


^see Appendix B for an introduction to the relevant aspects of this notation. 


















We see that the derivative operators on equations ([^ & Q provide m-raising 
and lowering operators on the radial basis. Equations 0 & ( |18[ ) in many ways 
resemble the same operators applied to Bessel functions, 


_d 

dr 


m 




dr 


m 




( 21 ) 


This is no accident. Both Bessel functions and the scaled Jacobi polynomials 
behave similarly near the origin. Alternatively, for many reasons Bessel func¬ 
tions are not as useful for representing the solutions to arbitrary PDEs [5]. In 
general, Bessel-function solutions achieve only algebraic convergence rates. For 
example, / = — 1 satisfies V^/ = 4, /(r = 1) = 0. This exact polynomial 

solution only decays as in a Bessel-function series expansion. This 

behaviour results from a mismatch between the higher-order derivatives of an 
exact solution and a given approximating expansion; e.g., the third derivative of 
— 1 vanishes identically, which does not happen for any finite Bessel function 
series. The same principle favours orthogonal polynomials over special functions 
in many other cases pK] . 


3. Tensor calculus 

We start with a scalar function 

OO 

f{r,9)= (22) 

m——oo 

Using the bra-ket notation, we expand the radial dependence of the Fourier 
components in terms of the weighted Jacobi polynomials 

OO 

fm{r) = (0,TO,r|/™) = ^/m,„(3°’™(r). (23) 

n—0 

We expand all scalars, vector, and tensor components in terms of the fc = 0, and 
appropriate m basis. As we differentiate a given component, the basis changes its 
m index, and the k index increments upward. A later section defines operators 
to convert between different fc, m bases. 
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The covariant derivative does not couple the azimuthal Fourier modes. We 
can therefore write the action of its radial and azimuthal components on an 
individual Fourier mode as 


V./™(r) = ^fmir), Vgfmir) 


Using equations ( [l7| ) & (18 1 , 


Vr/m(r) = ^(l,TO-l,r|£i \fm) + -^ 

= -^{ 1 , 171 - l,r\D-\ fm) - 


{l,m + l,r\D+\fm ), 
( 1 ,TO + l,r jiD+l fm )■ 


(24) 


(25) 

(26) 


Introducing a spinor basis simplifies the computations of vectors and tensors. 
The spinor basis elements are 


e± = ^ (Cr T ieg), 


(27) 


where er,eg represent the unit vectors in the r, d coordinate directions respec¬ 
tively. A 2 X 2 unitary matrix transforms between the components of a vector 
in the r, 9 basis and the ± basis. In the spinor basis, the covariant derivative 
becomes 


V — — e_V_ -1- 


(28) 


where acting on individual scalar functions, 




_d 

dr 


am 

r 


with a = ±1. From standard vector calculus, 


(29) 


V 7 - Cf — Vt" — 0, Cf — 6 ^, CQ — C-p. (^0) 

r r 

The spin basis diagonalises the connection coefficients such that 

(31) 

where cr,/r = ±1. From the standard three-dimensional vector algebra of 
Cj., eg, 63 , we can easily deduce 


e+X e_ = 163 , e+• e_ = 1, e+• e+ = e_ • e_ = 0, (32) 
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where the dot product represents simple component-wise contraction, *.e., not 
the proper complex-valued inner product. Further, 

63 X e± = zLie±. (33) 

We note that 63 denotes any locally orthogonal third unit vector. For example, 
it could represent the axial direction in cylindrical coordinates, or an additional 
angle in a torus. 


Equations (29) & (31) imply the following formula for the covariant derivative 


of an arbitrary tensor component 


Vcre^{k,m +n,r\ = /c, m-I-r I-I-fc, to-I-/ i, 1 

V2r 


d airn + /i) 


\/2 idr r 

= ei^(^k + l,m + + a,r. 


(k,m + fj,,r\ 


(34) 


Moreover, this formula holds if fi represents a multi-index such that 






s 

and = 

2=1 


(35) 


where fii = ±1 for each i. It fi G S = [—1,-|-1]* represents an s-component 
multi-index, then 


Tm(r) = '^efj,{k,m + n,r\TI^) (36) 

mgs 

represents a given rank-s tensor. The covariant derivative becomes the rank- 
(s -I- 1) tensor 


VT„(r) = ^ ^ ecre^( A:-I-1,TO-I-/r-I-cr,r Ill'll ). (37) 

(T — zt 1 fl^S 


In equation (37), k represents whichever basis the tensor began with. It need 
not equal the tensor rank, s. In most cases fc = 0. 


Computing the gradient of a scalar 

V/m(?') = e_( 1,TO - I,r |£)“|/m )-k e+( 1,TO-k l,r |l?+|/m ). (38) 


II 









We can further compute the second covariant derivative of a scalar 

VV/m(r) = e-e-{2,m-2,r\D~D~\fjn) + 
e+e-{2,m,r |D+£)“| fm) + 
e_e+( 2,m,r \D~D'^\ fm ) + 
e+e+{2,m + 2,r\D+D+\fm). (39) 

The operators commute in the sense that 

^k+l,m+l^k,m ~ ^ k+l,m-l^ k,m^ (4*^) 

where the k, m indices show that operators take on different meaning depending 
on their order. This implies the covariant derivative components commute; as 
they must for a flat space. Taking the trace of the second derivative gives the 
Laplacian of a scalar 

VVm(0 = 2{2,m,r\D+D-\fm). (41) 


4. Multiplication and Non-Constant Coefficients 


Many applications require multiplying scalar and vectors function by a non¬ 
constant, axisymmetric function of r. Here, we present a pair of operators 
representing multiplication without spoiling the bandedness of the matrix. In 
their discussion of vector formulations on disks, Sakai & Redekopp (2009) used 
vector components of the form rvm{f) |31j . These components behave like 
scalars. Two additional Jacobi polynomial recursions reconcile this formulation 


with our current analysis. Equation (11), and equations (12) & (14) together 


imply two distinct r-multiplication operators 

r(fc,TO,r| = (/c, m ± 1, r 


(42) 


where the dependency on k and m is dropped when inferred from context. 
R'^ has entries on the diagonal and first super-diagonal, while R~ contains 
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entries on the diagonal and first sub-diagonal. For m = 0, R = R~^. Equa¬ 


tions (A.6) & (A.9| show the individual elements associated with the operators 


in equation (42). 


Given the gradient vector in equation (38), 


f’V/m(r) = e-{l,m,r\R'^D \ fm) + e+{l,m,r\R D+\fm). (43) 


Both components of equation (43) contain the same m dependence as the scalar 
fmir). 


Considering equation (|^ and both ± parts of equation (42), 
z(^k,m,r\ = (^k,m,r \ {2R~R'^ — I), 


(44) 


which leaves /c,m unchanged. In the same style as equation (40), radial multi¬ 
plication commutes with itself, 




(45) 


where we note that the ± operators act on different m-bases depending on their 
order. 


The standard three-term recursion relation for Jacobi polynomials implies the 
symmetric three-term recursion for the radial functions 


Qt'-iir) + (r) = ^ Q^’-(r), 


(46) 


where 


A 


k,m 

n 


B 


k,m 

n 


m? — 

(2n -I- A; -I- m){2n + k + m + 2) 

2 n{n + k){n + m){n + k + m) 

2n + k + m\] (2n + k + m)^ — 1 


(47) 

(48) 


Therefore, 


Z = 2R-R+ - I 


(49) 
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represents the standard symmetric tridiagonal Jacobi matrix for Jacobi poly¬ 
nomials. The operation of multiplication by any purely radial function occurs 
via 


F{2r‘^ — l)(^k,m,r \ = (^k,m,r\F{Z). 


(50) 


We expand F(z) in terms of a truncated series of any polynomials, 


Njr 


F{Z) = Y^FnPniZ) 


(51) 


n—0 


An orthogonal polynomial basis allows the advantage of a matrix-valued three- 
term recursion of the form, 


Bn+lPn+l{Z) -I- AnPn{Z) + Pn-l{Z) — Z.Pn{Z) 


(52) 


Equation (52) need not represent any of the recursions in equation (46), but 


it can. For example, if we choose to expand in terms the same basis as the 
axisymmetric k = m = 0 modes, then 


An — 0 , Bn = 


(53) 


V Alp — 1 

But this is not an essential restriction. We could use Chebyshev polynomials 


for Pn{Z) if the series in equation (511 converges more rapidly in that basis. In 


that case, A„ = 0, = 1/2 for n > 0, and Bq = '^/y/2. 


Reducing the truncation degree, Np^ reduces the bandwidth of the multipli¬ 
cation operator matrices and can save signihcant computation time even for 
modest differences in the truncation degree. Z is a, tridiagonal matrix, therefore 
the bandwidth of F{Z) is 2Nf + 1, which is typically much smaller than the 
required size of the matrix. 


Given any choice in equation (52) and taking P_ 
algorithm to compute F{Z) directly via 


= 0, we can use Clenshaw’s 


Knf + 1 = Knp+2 = 0 


(54) 
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Kn — Fnl + 


1 


(55) 


-- {Z - AnI).Kn+i - 

-On+1 -On+2 


F[Z) = PoiZ)Ko, 


(56) 


where / represents the identity matrix, and Po{Z) = PqI is a constant times the 
identity matrix. Each iteration of Clenshaw’s algorithm increases the bandwidth 
by one, via multiplication by the tridiagonal matrix Z. If the coefficients vary 
strongly with latitude, then the multplication matrices are dense. But this 
happens in a gradual way. Most of the time, the coefficients do not vary nearly 
as rapidly as the solution itself. 


Multiplication by non-azimuthally symmetric functions must obey more com¬ 
plicated selection rules for the different m indices. This would create systems 
that couple both the 9 and r directions. We do not consider linear operators 
with m > 0 dependence. 


5. Basis conversion 


Equation (41) has a problem: the Laplacian of a scalar with k = 0 results in a 


scalar with k = 2. We must convert between the k bases in order to compare 
the input and output. 


Equation (121 implies a 2-band conversion operator such that 
(fc,m,r| = {k + l,m,r\Ck,m, 


(57) 


where the dependence on k and m is again dropped when inferred from context. 
C has entries on the diagonal and first super-diagonal. Like the operators. 


we omit the k^m indices when the context makes it clear. Equation (A. 121 


shows the individual elements associated with the operators in equation (57). 


Raising the k index allows converting to the basis that most naturally represents 
derivatives. This is the same principle that underlies sparse representations with 
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Chebyshev polynomials, i.e., 

^T^{z) = nUn-i{z), 2r„(z) = Un{z) - Un-2{z). (58) 

dz 

Many papers reinventing Chebyshev numerical schemes rely on this simple fact 
[571 m [H [TH [51 [in m Hn inn mm . chebyshev polynomials are special cases 
of Jacobi polynomials with a = 5 = —1/2 and a = h = 1/2 representing the 
first and second kind, respectively. Using greater flexibility in the a, b indices 
allows for an increased ability to construct sparse operators for a particular 
application. 

5.1. Heisenberg Algebra 

Radial multiplication, differentiation, and basis conversion satisfy simple alge¬ 
braic compatibility conditions. Differentiation and multiplication satisfy stan¬ 
dard commutation relations, but with the added element of conversion. 

The commutation relations satisfy the two-dimensional Heisenberg Lie algebra, 
() 2 . That is, 

[D+,R-] = [D-,R+] = V2C. (59) 

Keeping track of the different fc, m indices, 

^k,m^l^^,Tn ~ ^^+l,Tn±l^k,m ~ ( 60 ) 

Otherwise all other ± operators commute with each other and with C. Fig¬ 
ure shows the structure of the different elements of ()2 acting on the ( fc, m, r | 
state. The diagram omits the equivalent operators acting on the other nodes in 
the infinite {k,m) lattice. Traditional representation theory considers the Lie 
algebra of operators acting on a single Hilbert space. Here it appears that we 
are considering the unconventional situation of a Lie algebra of operators acting 
between Hilbert spaces; each with different k,m. Considering a total Hilbert 
space comprising the direct sum of all individual fc, m spaces rectifies this issue. 


16 



(fc + 1, m + 1) 


{k, m + 1) 


Figure 1: Diagram of operator structure for the (fc, m) indices. Every other node in the graph 
contains its own set of operators mapping to new nodes in a similar fashion. Only the right 
half of the diagram exists for m = 0. 


The operators also satisfy conditions apart from commutation relations. Taking 
the 9 component of equation (|43|) implies 

(61) 


1 


mC =-^ {R+D- - R-D+) . 
Taking the r component of equation ( |43[ ) implies, 

-ic=^{R-^D-+R-D^). 


V2 


(62) 


Equation (621 is equivalent to equation (15) from Matsushima & Marcus (1995) 

Hi. 


The Heisenberg algebra allows the automatic reformulation of equations at the 


numerical level. Equation (59) gives a simple and reliable recipe for commuting 


radial non-constant coefhcients past derivatives using the numerical operators. 


5.2. Boundary restriction 


Solving PDEs on the unit disk usually requires appropriate boundary conditions 
at r = 1. No boundary at r = 0 is needed, nor allowed. To impose the boundary 
condition at r = 1 we use 


g*’'"(r = l) = W2(2n + m + fc+l) 


k 


n + m + k 

k 


(63) 
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which follows from equations (131 & (15). We can build higher-order boundary 


restriction operators out of equation (63) and various differentiation formulae. 
For example, 




2n{n -I- fc -I- 1) -I- m{2n -|- fc -|- 1) ^ 


It is convenient to modify the basis so that the boundary condition only ap¬ 
pears in the first element. This leads to banded systems, as opposed to almost 
banded systems, which may have dense rows corresponding to the boundary 
conditions. Many black-box linear algebra packages can solve banded systems 
very efficiently, while almost-banded systems require specially designed schemes, 
e.g., the adaptive QR method [53], or use of the Woodbury formula. 


We describe how to modify the the basis to accommodate Dirichlet boundary 


conditions, corresponding to (13). The adjoint of the conversion operators gives 
a fc-lowering operator. 

(1 — r^)( fc, TO, r I = ( fc — 1, TO, r (65) 

where has entries on the diagonal and first sub-diagonal. This follows from 


a straightforward application of equations (11) & (14). Equation (A. 15) gives 


the explicit indices of equation (65). The fc-lowering operator satisfies the com¬ 
mutation relations, 


[C\D^] = V2R^. 


( 66 ) 


Otherwise commutes with all other relevant operators. The operators 
are the adjoints of each other. Adjoints of the L>^ operators exist, but provide 
little apparent use for the types of problems we consider here. 


Define 


B = + a ^|0)(0|, where a=(fc,TO,r = l|0), (67) 

where ( 0 | represents the unit basis element corresponding to the n = 0 mode. 


Equation (67) defines a two-band lower-diagonal operator. 


18 













As before, 


fm{r) = {k,m,r\fm) (68) 

If we change variables such that, 

|/„) = B\g^) ( 69 ) 

then 

fm{r) = {k,m,r\C^\gm)+ )- -rnTY ( ^ 5m ) 

' ' ^ /c, m, 1 I 0 y 

= (l-r^)(fc + l,m,r|g^)+ / f ’ j n ( 0 | 5m )■ (70) 

' ^ /c, m, 1 I 0 y 

Therefore 

fm{r = l) = (0|g™). (71) 

The B operator localises the boundary-restriction operator to the first mode of g. 
This contrasts with a dense row operation for the original variable. This comes 
at the expense of a single band in the bulk equation operators, but produces an 
entirely banded system. 


For the derivative, because for every individual element ( fc, m, r = 1 | >0, there 
exists a diagonal operator A such that 

fmir = l) = (fc,m, l|A|/™) (72) 

This also extends to more general left-hand side restriction operators. In the 


case of the derivative, equation (64) implies the diagonal elements 


( n I A| n ) = 


2n{n -I- fc -I-1) -I- TO(2n -|- fc -|- 1) 
k + \ 


Therefore, if 


then 


fm) = A ^B\g^), 


(73) 


(74) 


/m(’’= 1) = {k,m,l\B\g,ri) = (0|g„). 


(75) 
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The same principle applies for any A, not only for the derivative. Julien & Wat¬ 
son (2009) demonstrate many of the advantages of a sparse transformation to 


a Galerkin basis for general boundary operators M- Equation ([73| shows that 


A„=o = 0 for TO = 0, otherwise A„ > 0. This precludes forming equation (741 
for n = TO = 0. Neumann boundary conditions fail to determine this particular 
mode, which requires additional information. In many applications, setting this 
mode requires fixing a physically irrelevant (but nevertheless mathematically 
essential) gauge condition. 


5.3. Quadrature 


All the above analysis considers field variables in the spectral coefficient do¬ 
main. This allows for efficient differentiation operations. Given a set of spectral 
coefficients, we need the option to construct the function on a radial grid. We 
also need to compute spectral coefficients from data on a grid. In general, we 
need each to mode on the same radial grid. This allows constructing solutions 
on the full r, 9 disk. We consider an expansion in terms of a finite number of 
radial modes, 

Nk.^ 

fm{r) = ^ Q^’'"(r)(fc,TO,n|/m), (76) 

n—0 

where Nk,m depends on both k Sz m. We use Gauss-Legendre quadrature on 
the grid Zi = 2r’f — 1 such that 

Nr- 

{k,m,n\fm) = '^Qt''^{r^){l - r‘f)'^Wif^{r,) (77) 

i=l 

for some W- Gauss quadrature is exact for all z polynomials less than degree 
2Nr, or r polynomials less that degree 4AV. In the r variable, 

deg [Q^’'"(r-)(I - r")Vm(r)] = 2(n + to + fc + Nk,m) < 4iV, (78) 


where n < Nk^m- Therefore, 


fc,m 


Nr — 1 — floor 


k + m 
2 


(79) 
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yields exact quadrature on a Gauss-Legendre grid. Gauss-Legendre quadrature 
is based on the fc = 0, m = 0 case. For non-integer starting k values, Wi, ri would 
result from a more general Gauss-Jacobi quadrature, but the basic principle 
remains the same. In practice, we never need to compute equation ([77]) with k > 
0. Instead, we can project any function on to the k = 0 spectral coefficients 
and convert the result up the desired k using the C conversion operator. 


The absence of a fast transform is the largest downside to the Jacobi-based 
method compared to a Chebyshev scheme. In extreme cases of large Nr, where 
a fast transform would show signihcant gains, Chebyshev schemes struggle to 
represent the severe analytic behaviour of the solution near r = 0. For example, 
if m is an even integer, with p = mj^l, and k = 0, then 


- 1 ) = 

"+P n+p , 2z \ri\rn+p\/i+n+p\ . , w+„+p 

E Lli* ■ " - I. - ( 80 ) 

j=0 i=max{j,p} 




where T 2 j{r) represents a Chebyshev polynomial of even degree 2j. The primed- 


summation notation in equation (801 represents multiplying the J = 0 term of 


the series by 1/2. The series in equation (80) contains Chebyshev modes up to 
degree 2n-|-TO. A similar relation exists for odd m values. Given a relative error 
tolerance e, in some cases these series coefficient will become smaller than the 
allowed error before the end of the exact series. However, even this does not gain 
much advantage. In the regime m ~ n, a Chebyshev series expansion of Q^^(r) 
requires nearly all of the 2n+m modes to achieve an accurate representation. We 
note that these considerations ignore the possibility of a low-rank approximation 
to specific functions, which can also exploit a fast transform; see [3^ . 


Also, because Nk^m decreases with larger to, a Jacobi-based scheme requires 
fewer total modes to obtain uniform resolution over the entire disk. That is, 


after summing equation (79) over |to| < N 0/2 we find 


Ntot. 



(81) 
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compared to simply « NrNg for a Chebyshev radial scheme. The savings are 
maximised for Ng « diV^, with A^tot. ~ NrNg/2. 

Two final mitigating factors arise from (*) the somewhat-large cost of fast trans¬ 
forms for moderate sizes, and (ii) the easy parallelisation of matrix transforms 
(i.e., quadrature). Tablecompares the timing of a fast Chebyshev transform 
(FCT, based on the type-II fast cosine transform), versus quadrature computed 
with a matrix-vector multiply (MMT). The tests show that while the FCT gives 
better asymptotic scaling, the MMT does better at moderate sizes. The speed 
crossover occurs somewhere between 64 and 128 modes; a reasonable size for 
many applications. The timing difference remains less than an order of magni¬ 
tude all the way up to 1024 modes. Lastly, a quadrature-based transform uses 
the orthogonality of the eigenfunctions to compute the inverse transform from 
the transpose of the forward transform. That is, both the forward and inverse 
transform represent ^ operations. However, because of the matrix-vector 

structure, the operations actually entail Nr totally independent vector-vector 
dot products. This parallelises quite easily. 


Nr 

: 32 

64 

128 

256 

512 

1024 

2048 

FCT (/Ltsec) 

: 0.94 

1.97 

4.42 

8.56 

17.8 

40.7 

81.3 

MMT (/Ltsec) 

: 0.43 

1.75 

9.19 

30.1 

117. 

462. 

1905. 

MMT/FCT 

: 0.45 

0.89 

2.08 

3.51 

6.59 

11.3 

23.4 


Table 1: Timings for Fast Chebyshev/Cosine Transform (FCT) versus a Gauss-Legendre 
Matrix-Multiply Transform for different resolutions. The timings use SciPy (FCT) and 
NumPy (MMT) on a Late-2013 Mac Pro 2.7 GHz 12-Core Intel Xeon E5. Each timing 
represents the average over 30,000 forward and inverse transforms of random data; only the 
actual transforms factor into the average. 


6. Examples 

This section considers a series of practical examples. Example 1 demonstrates 
the operator methods on a simple scalar equation. Example 2 demonstrates the 
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same for a vector set of equations. Example 3 concerns non-constant coefficients 
in a vector system. Example 4 considers a problem with full 0, r-dependence 
(z.e., with a large number of m values) and considers computational perfor¬ 
mance. 

The first three examples entail finding the generalised eigenvalues and eigenvec¬ 
tors for an analytically tractable, and/or previously published numerical test. 
We compute eigenvalues because of the demanding nature of these problems. 
Eigenvalues and eigenvectors give a significant fraction of all the useful informa¬ 
tion about a given system. For example, the eigenvalue problem determines the 
stability of general time-stepping schemes. In some cases, defective numerical 
schemes can produce large spurious eigenvalues in the wrong part of the com¬ 
plex plane [Samuil]- Such methods can give accurate solutions for a single 
linear solve with specific right-hand side terms. But rogue eigenvalues in linear 
operators can render a time stepper useless after a small number of iterations. 
A small random projection onto a bad eigenvector will diverge exponentially 
upon iteration. Also, eigenvalues in the wrong part of the complex plane can 
prevent otherwise conserved quantities from remaining constant. 

In each case we build the operator matrices by inspection of the original equation 
and use standard linear algebra libraries for the actual calculations. Apart 
from example 4, we code all the examples in Python, and take advantage of 
the features within NumPy and/or SciPy. We produce Jupyter notebooks, and 
scripts for each case and make them available online within the Dedalus projecij^ 

Example 4 is coded in the ApproxFun.jl Julia package [32 [33], to utilize its 
implementation of the adaptive QR method. 


^dedalus-project.org 
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Example 1 (Bessel function eigenvalue problem):. We solve the following eigen¬ 
value problem 


l_d 

r dr 


r^Jminr) 

dr 


9 ^ 


(82) 


with the Dirichlet boundary condition, 


dm (^) — 0- 


(83) 


Many available software packages compute Bessel functions and their corre¬ 
sponding zeros. We therefore use this example as a numerical test. 


We solve the generalised eigenvalue problem, 

2D-D+\fm) = -K^CC\U), (84) 

where CC denotes applying two (different) conversion operators. The left-hand 
side operator only contains a single off-diagonal band. Each C operator only 
contains a single off-diagonal band. The product of the two is tridiagonal. This 
contrasts with the penta-diagonal system for this problem in |18) . 


The boundary restriction operator in equation (631 replaces the last matrix row 
on the left-hand side. A row of zeros replace the last matrix row on the right- 
hand side. After constructing the appropriate right- and left-hand side matrices, 
we use the NumPy eigenvalue solver to find the eigenvalues and eigenvectors 
in coefficient space. We sort the results from lowest-to-highest in eigenvalue, 
and transform the coefficient-space vector to a radial grid; see §5.3 for details 
regarding the spectral-to-spatial transforms. 


Figure shows the relative error in the eigenvalues for an to = 50 solution set 
with Nr = 500 radial modes. The plot is typical of the behaviour for other to 
values. We choose this case to examine closely because it gives a good illustra¬ 
tion of the ability of the method to capture the large dynamic range near r = 0. 
However the method’s performance near the origin is independent of the to 
value. The behaviour of the numerical eigenvalues is typical of orthogonal poly¬ 
nomial solutions of boundary-value problems, including those using Chebyshev 
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polynomials. Typically, given a value of Nr, approximately the first 1/2 to 2/3 
of the spectrum will contain values with virtually no error; apart from machine 
floating-point error, and a little from the eigenvalue solver routine. In figure 
this region spans from mode 0 to about mode 300 out of 500 modes. After this 
point, the exact solutions begin to outstrip the inherent finite-degree polynomial 
approximation. This error increases rapidly, with the largest eigenvalues scaling 
such that 

|Ak„| ~ (85) 

Like Chebyshev schemes, this behaviour is typical for orthogonal polynomial 
solutions of two-point boundary eigenvalue value problems. Figure implies 
that obtaining accurate solutions to Poisson’s equation, 

VVm(r) = s™(r) ^ 2D-D+\fm) = CC\sra) (86) 

requires that Smir) does not project onto the inaccurate part of the spectrum 
for a given resolution. This implies that the resolution needed to compute fmif) 
is roughly 3/2 to 2 times larger than needed to resolve Sm{r)- But after reaching 
the critical number of modes, the error is effectively machine precision. This 
same principle holds for most types of equations. 

Figure shows the eigenfunction behaviour for the n = 200 eigenvector for m = 
50. This corresponds to K 200 ~ 707.447066905. The highly oscillatory solution is 
accurate to roughly 13 digits mostly uniformly throughout the domain. Figure]^ 
shows the spectral coefficients for the same eigenfunction. The magnitude of the 
coefficients decrease dramatically after approximately 300 modes. This accords 
with the similar behaviour in the the eigenvalues. Trying to represent this mode 
with fewer than 300 to 400 modes would incur significant error. Beyond this 
number, more resolution contributes little to the fidelity of the solution. 

Finally, we note that the matrices for radial component in polar coordinates 
are sparser than the Cartesian Chebyshev analogue in one dimension. Cheby¬ 
shev discretisation applied in polar coordinates must account for the geometric 
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eigenmode number, n 

Figure 2: Plots of relative eigenvalue error from examples 1 Sz 2. Each plot compares the 
numerical eigenvalue solution with the analytically computed value for 500 individual sorted 
modes. Panel (a) corresponds to Bessel function example 1 for m = 50. Panel (b) corresponds 
to the inertial waves example 2 m = 1. Each case contains roughly 70% ’’good” values with 
high accuracy, and 30% ’’bad” values with significant error. In both cases, the bad values still 
lie in the same region on the complex plane as the true solutions. The bad values correspond 
to small-scale modes that outstrip the fixed resolution. In the Bessel function example, small- 
scale modes correspond to large k. Hence the error become large for large mode numbers. 
For inertial waves, small-scale modes correspond to slow frequencies with |cj| 5^ 0. The more 
accurate large-scale modes correspond to |a;| 5^ 1. Therefore, the significant errors on panel 
(b) occur for intermediate mode numbers; see Table § 
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0 

-2 



Figure 3: Plots for the 200th eigenfunction from Example 1 for m = 50, and Nr = 500. 
Panel (a) shows a log-log plot of radial dependence near the origin. The solid blue line shows 
the computed eigenfunction. The eigenfunction should scale as ~ The dashed black 

line denotes this scaling. Panel (b) shows the eigenfunction across the entire radius. The 
solution vanishes at r = 1. Panel (c) shows the error between the computed solution and an 
analytically computed Bessel function. The characteristic error is Pi 2 x 10“^^, and mostly 
uniform over the domain. 
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Figure 4: Plot of the spectral coefficient amplitude of the 200th eigenfunction from Example 1 
for m = 50. In general, the nth eigenfunction requires P=:2n spectral coefficients to resolve. 


terms with non-constant coefficients. Using operators that conform to the disk 
geometry renders systems such as equation (82) effectively constant coefficient 
problems. 


Example 2 (Inviscid inertial waves in a upright rotating cylinder):. The fol¬ 
lowing example highlights vector-scalar coupled system with a non-trivial an¬ 
alytical solution. It therefore makes a good system for numerical comparison. 
We consider 


iu Vm{r) - 1 - 63 X Vm{r) + '^Pm{r) 

= 0, 

(87) 

iujV ■ Vm{r) + a^Pmir) 

= 0. 

(88) 


The two-dimensional velocity vector and pressure, Vm{r) = {r)e+ + v^{r)e- 

and Pmif) represent a full eigenstate. The parameter uj represents a frequency 
eigenvalue. The parameter a represents a mode aspect ratio in the third di¬ 


mension. We solve equations (87) & (88) with the normal-velocity boundary 
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n : 

0 

100 

200 

300 

400 

500 

Kn {m = 50) : 

5.71el 

3.91e2 

7.07e2 

1.02e3 

1.33e3 

1.64e3 

Wn (to = 1) : 

-2.14e-l 

-3.17e-3 

-1.23e-3 

1.20e-3 

3.12e-3 

3.19e-l 


Table 2: Approximate analytically computed eigenvalues versus mode number. The approx¬ 
imate values are intended to give scale to the bottom axes on Figure]^ txn increases mono- 
tonically from large-scale modes to small scale. For ojn, small-scale modes correspond to 
the middle of the sorted eigenvalue spectrum, and large-scale modes correspond to the val¬ 
ues closer to unity. Inertial waves break reflection symmetry. Hence the lack of symmetry 
u —UJ. 

condition 


Greenspan (1968) 


G • v^{r = 1) = 0. (89) 

discusses at length the derivation and physical interpre¬ 


tation of equations (87) & ( [88| . We take the system as given. Inertial waves 
make a good test problem because of the non-trivial frequency constraint 

- 1 < w < 1. (90) 


As with the Bessel function problem, numerical eigenvalue solvers often produce 
a subset of values that differ significantly from their true solution. It is, however, 
a sign of a good scheme if all values are maintained in the same region of the 
complex plane as the true result. As discussed above, a poorly designed scheme 
can produce very large spurious eigenvalues with a different complex character. 


An analytical solution to equations (87)-(891 results from collapsing the vector 
system to the scalar problem 


V72 / ^ (l-w2)a2 

V Pm{r) = —-2- -Pm{r) 


(91) 


with 


d m 

W 3 -Pm(r)-I - Pm[r) = 0 at r = 1. 

dr r 


(92) 


Equation (91) requires a solution in terms of Bessel functions, Pniir) = Jminr). 
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Equation (921 requires 


= 


(1- w 2)^2 


KU!J!^{k) + mJm{K) = 0 . 


(93) 


We solve equation (93) with a standard numerical root finder. 


Recall equation (33), 63 x e± = ±ie±. This implies the following block eigen- 
system 


(94) 


c 

0 

D+ 


c 

0 

0 

0 

-c 

D- 

1 ^ 

0 

c 

0 

0 

0 

a^C 


D- 

D+ 

0 


where the eigenstate is the vectorisation of the velocity and pressure, 

T 


) = l^m) V'Vm) \Pm) 


(95) 


Note that absorbing a factor of v — 1 into the velocity field renders equation (94) 
purely real. The boundary conditions are 


(0,TO-fl,l| (0,m-l,l| 0 |S) = 0 


(96) 


This boundary-restriction operator replaces the last matrix row on top blocks 


of the left-hand side of equation (94); the blocks corresponding to u'*'. A row 


of zeros replace the last matrix row of the top blocks on the right-hand side. 


We solve the numerical discretisation of equation (94) with the same eigenvalue 
software as example 1 . 


Like example 1, the inertial waves problem only allows one boundary condition 
for a system with (apparently) two radial derivatives. This results from the 
coordinate singularity causing a reduction in the differential order at r = 0 . 
This manifests in the matrix system because D~ is a non-singular diagonal 
matrix, and hence requires no additional conditions to invert. Alternatively, 
the D'^ matrix contains a null space, and requires a boundary condition. We 
insert the radial-velocity boundary condition in the block corresponding to the 
I ) variable. 


30 

















Figure shows the relative eigenvalue error for the inertial wave problem with 
TO = 1, a = 1. Like the Bessel function example, there is a significant region 
with very high accuracy, and a region of the spectrum with significant error. 
Both of these regions occupy roughly the same proportions as the respective 
regions in the Bessel function problem. The major distinction is that the error 
occurs in the middle of the spectrum, rather than on one end. This results from 
the fact that small-horizontal-scale inertial wave oscillate slower than domain¬ 
filling modes. The inaccurate part of the spectrum results from under resolving 
small-scale modes in all cases. In both example 1 & 2, all the eigenvalues exist 
in the correct part of the complex plane; on the positive real axis for example 1, 
and on the real axis between ±1 for example 2. While it is not possible to avoid 
inaccurate values for some modes, remaining on the real axis is not a guaranteed 
feature of a numerical approximation. Our scheme succeeds well in this respect. 

We discuss to = 1 modes in this second example, but note that all modes 
preform well, just like in example 1. In the higher-m cases the spectrum becomes 
squeezed to a smaller region between ±1, otherwise there is little difference. The 
eigenfunction error corresponds to the error in the eigenvalues. 

Example 3 (Linear perturbations of Hagen-Poiseuille pipe flow):. Hagen and 
Poiseuille first studied pipe flow experimentally in the 1830s [53] . This problem 
has guided studies investigating the onset of turbulence since the pioneering ex¬ 
periments of Reynolds in the 1880s [5S]. A lot of work continues today regarding 
the nonlinear dynamics of this system [27j . After much investigation, a large 
sector of the turbulence-modelling community believes all linear perturbations 
in Poiseuille pipe flow eventually decay exponentially in time. Finite-amplitude 
perturbations likely produce the rich turbulent behaviour observed in labora¬ 
tory experiments. This problem interests us because it poses a difficult numer¬ 
ical challenge; independent of the physical significance of the linear dynamics. 
Meseguer & Trefethen (2003) |5D| discuss the linear problem at length. The 
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eigenvalue problem for linear perturbations is: 


Aum(r) -1- Vpm(r) -|- Lvm(r) 

= 0 

(97) 

Vm{r) + iapm{r) + Lw^{r) 

= 0 

(98) 

V • Vm{r) +iawm{r) 

= 0. 

(99) 


where 


L = iaW{r) — . 


( 100 ) 


For the background flow in the stream-wise direction, we choose the pressure- 
driven parabolic profile 


W{r) = 1 — with W'{r) = —2r. 


( 101 ) 


The held Um(r) represents the two-dimensional vector velocity in the pipe cross 
section. The scalar helds Wmir), and Pmir) represent the velocity in the stream- 


wise direction, and pressure respectively. The Laplacian implicit in equation (1001 


represents the scalar Laplacian in equation (981, or the (more complicated) vec¬ 


tor operator with curvature terms in equation (97). The advantage of the spinor 


approach is that these present the same level of complication in the numerical 
implementation. The parameters a, A represent the mode wavenumber along 
the pipe axis, and the complex-valued growth rate respectively. The viscos¬ 
ity parameter v = 1/Re, where Re represents the how Reynolds number. We 
require no-slip boundary conditions on the outer wall 


,(r = 1) = Wm{r = 1) = 0. 


The full eigenstate now contains four held components 


) - km) km) km) km) 


n T 


( 102 ) 


(103) 


We solve the generalised eigenvalue problem for the complex-valued growth rate, 
A, 




( 104 ) 
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where 


L = 


7 ^ = 


L 

0 

-V2CCR- 

D- 

CC 0 
0 CC 
0 0 
0 0 


0 

L 

- -V2CCR+ 
D+ 

0 0 
0 0 
CC 0 
0 0 


0 

0 

L 

iaC 


CD+ 

CD- 

iaCC 

0 


(105) 


(106) 


and L = iaCC{I — R'^R ) — 2vD'^D + va^CC. The boundary conditions 
become 


) = 0. 


(107) 


(0,m + l,l| 0 0 0 

0 (0 ,m-l,l| 0 0 

0 0 (0,m, 1| 0 

As with all other examples, we note that the operators in the different blocks of 


equations (105) & (1061 acquire different (fc, m) values depending on the field 
components. 


In addition to the eigenvalue problem, we solve for the parallel-component of 
the vorticity, wa, and stream function, '(/). 

= V • (ea X Um) = ■i(V_u+ - V+u“) = wa, (108) 

which translates to 


2D- D+\ijm) = iCD-\v^) - iCD+\v;;,). 


(109) 


Often for an incompressible flow, the parallel velocity and stream function pro¬ 
vide the dynamical variables in a higher-order-derivative scalar formulation of 
the original vector system. This example requires no special treatment to en¬ 


force the incompressibility constraint in equation (99), which stands on equal 
footing with the other equations in the system. 
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Figure shows full-disk eigenmodes for the first- and second- slowest decaying 
modes for Re = 10^/^, a = 1, and to = 1. The plots show the four primitive 
dynamical variables p,Vr,V 0 ,w along with the two ex post facto derived fields 
uj 3 ,ip. Below each full-disk image {x^ + < 1), we show a zoomed-in square 

image on the region max(|a;|, |y|) < 0.2. The zoom-in plots show no signs of 
numerical singularities for at the origin. The scalar fields, p,w,u! 3 ,tjj remain 
smooth, and the vector quantities show the correct parity. Figure shows ad¬ 
ditional full-disk eigenmodes for to = 5, and to = 12. These plots show that 
the higher-TO still behave as required near the origin. Traditional methods {e.g., 
Chebyshev bases) can cope with 1 ow-to values, but fail for to 3 > 1. Nonlin¬ 
ear simulations at high Reynolds number require a large range of azimuthal 
wavenumber to resolve the resulting turbulent states. For these situations, a 
scheme with a correct accounting for the origin is essential. 

Table shows a list of computed growth rates for different to values for a = 
1, and two different Reynolds numbers. In each case, we report the fastest 
growing (equivalently slowest decaying) modes of the centre- and wall-localised 
branches. The to = 1 values for both Reynolds numbers match the values 
reported in Meseguer & Trefethen (2003) [5D] to the previously reported 11 
digits. We tested all other values reported in m- We match all cases. We 
report additional (higher to) values for potential future numerical comparison. 
The individual modes reported in table correspond to the same eigenmodes 
shown in figures [^ & |^ For most cases, the first centre mode decays slower 
than the first wall mode. But this is not always the case. These modes can lie 
immediately adjacent to each other in the spectrum, but not necessarily so. For 
Re = 10^, we converge to the reported solution after « 50 radial basis elements. 
For Re = 10^ we converge after « 200 radial basis elements. In both cases the 
convergences is independent of to. 
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Figure 5: Full-disk images for the first two modes for pipe-flow with m = 1 and Re = 3162, 
and Nr = 200 radial points. Panel (a) & (b) shows the first- and second-slowest-decaying 
mode respectively. In each panel, the square plots show zoom-in plots for |a;|,|j/| < 0.2. The 
behaviour at the origin remains regular for each field. 
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Figure 6: Additional full-disk images for different pipe-flow modes with Re = 3162, and 
Nr = 200 radial points, (a) corresponds to m = 5, n = 1, (b) corresponds to m = 5, n = 2, (c) 
corresponds to m = 12, n = 1, (d) corresponds to m = 12, n = 5. Zoom-in plot are omitted, 
but the behaviour at the origin remains smooth for each field. 
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m 

type 

n 

Re 

Real(A) 

Imag(A) 

1 

centre 

0 

10 ^ 

-0.0227049145535 

0.951481194735 

1 

wall 

1 

10 ^ 

-0.0472321995947 

0.273788709331 

5 

centre 

0 

10 ^ 

-0.0725274157946 

0.898561158159 

5 

wall 

1 

lO'i 

-0.0793504734563 

0.247410847332 

12 

wall 

0 

10 ^ 

-0.0948648867252 

0.144951983763 

12 

centre 

3 

10 ^ 

-0.170456145014 

0.800901547889 

1 

centre 

0 

10 " 

-0.000721091206991 

0.998464685977 

1 

wall 

15 

10 " 

-0.00748956875998 

0.0303389812102 

5 

centre 

0 

10 " 

-0.00229096203822 

0.996790918537 

5 

wall 

15 

10 " 

-0.00855398926555 

0.0148836399355 

12 

centre 

0 

10 " 

-0.00538731680888 

0.993703412087 

12 

wall 

5 

10 " 

-0.00784725003139 

0.0296167267785 


Table 3: Pipe-flow eigenvalues computed in example 3. In each case, the “centre” and “wall” 
modes represent the slowest-decaying solution with eigenfunctions localised near the origin or 
near the boundary, respectively. The number n represents where a given mode occurs in a 
sorted list of growth rates. 

Example 4 (Forced Helmholtz equation):. We consider the forced Helmholtz 
equation 

(V^ + k2)/ = s and f{r = l,e) = g{e) (110) 

with s{x,y) = and g{x,y) = ycos{10x). When k is very 

large, the solutions become increasingly oscillatory; see Figure The zoom- 
in plot on the right demonstrates that no artihcial singularity arises at the 
origin, despite the complicated nature of the solution. We choose this problem 
to demonstrate the performance of the methods on a challenging and natural 
problem requiring simultaneous r, 9 dependence. 

The bandedness of our representation implies the complexity of the solution is 
0{MN) operations, where M is the number of 9 modes needed to resolve g{9) 
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Figure 7: Contour plot of the calculated solution to the forced Helmholtz equation for k = 60, 
with s(x,y) = £—(^—0-4,) —{a—0.3) g(x,y) = j;cos(10a:). No artificial singularity appears 

near the origin (marked with small dot in the zoomed-in image). 


and s(r, 9) in the 9 direction, and N is the optimal number of polynomial coeffi¬ 
cients needed to solve to a given tolerance in radius. In Figure|^we demonstrate 
the linear complexity for two different resolved right-hand sides. We approxi¬ 
mate, -(y-O-3) 25 9 modes and 11 radial polynomial coefficients 

per mode, while cos(40 sin(3a;2/) -I- 1) requires 149 9 modes and 60 polynomial 
coefficients per mode. We use the adaptive QR method implemented in the 
ApproxFun.jl Julia package [531 [21]. This adaptively determines the optimal 
number, N, of polynomial coefficients needed to resolve the solution for each 


9 mode. Equation (110) presents a challenge because the solution can become 


much more structured than the forcing for k ^ 1. At the high rage of k, we 
solve an equation with over 25 million unknowns in less than 25 second^ 


We also solve a screened Poisson equation similar in form the equation (1101 only 
with < 0. Rather than small-scale oscillations, the solutions tend to form 


^Timings on a 2012 MacBook Pro, with an 2.7 GHz Intel Core i7. 
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(b) 10' 




Figure 8: (a) Time taken to calculate solutions to the Helmholtz equation for two different 
forcing terms (green versus blue) with g{x^y) = :?/cos(10ic). (b) Total number of coefficients 
needed to resolve the solutions. The asymptotic slope becomes linear in k after an initial fixed 
costs needed to resolve the right-hand side forcing, as indicated by the background dashed 
grid. 


sharp boundary layers near r = 1. The performance efficiency and accuracy of 
the methods treats this case equally as well as the Helmholtz problem. 


7. Conclusion 

This paper produces a new set of algorithms for numerical computations on 
the unit disk. Choosing a hierarchy of bases built out of Jacobi polynomials 
allows for solving a range of important differential equations. We show a range 
of time-independent examples. Extending to time-evolving system only requires 
conducting the type of calculations presented here in the iterative fashion. 

The methods described here explicitly construct basis elements that respect 
the coordinate singularity in polar coordinates at r = 0. Beyond the disk, the 
type of calculus presented here will exist for other geometries with coordinate 
singularities. The two-dimensional disk represents only a prototypical example. 


39 











































The Laguerre class of polynomials allows representing functions in polar coor¬ 
dinates for the whole two-dimensional plane. In particular, 

lim = (Ill) 

i ?—>00 

Many of the formulae present in this paper extend to the infinite plane for 
fc —> 00 with only minor adjustments for the Gaussian factor. In this limit, the 
conversion matrix becomes the identity, C —>■ /. 

Two-dimensional spheres and solid balls also display similar behaviour as the 
unit disk. The sphere effectively contains a local disk at both poles. Spherical 
harmonic functions naturally resolve the coordinate singularities at the poles 
for scalar functions. However, a wider class of Jacobi polynomials allows the 
use of sparse vector calculus for higher-order tensor helds. In addition to the 
poles, the centre of a solid sphere requires relations completely analogous to 
equations Q & ([^ only for the spherical harmonic degree [55] . In an upcoming 
publication, we show how to organise sparse tensor calculus in three-dimensional 
spheres. This includes algebraic structures similar to those found here for the 
disk. 

The algebraic-analytic link is an important theme this current paper. As a 
general principle, the singularities of a PDE, either coordinate-based or phys¬ 
ical, dictate what types of bases are needed to represent meaningful solutions. 
When followed carefully, this guiding principle should allow the construction of 
numerical methods more efficient and straightforward to implement than more 
traditional schemes. Spherical geometry exhausts the parameter variation found 
in Jacobi polynomials. However, many more lesser-explored sets of orthogonal 
polynomials exist, with a large variety of symmetry properties. This includes 
discrete sets. We encourage the future development of numerical methods based 
on the guiding algebraic-geometric properties of a given general problem class. 
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A. Matrix entries 

In this appendix Sij represents the Kronecker delta. We show all matrix coef¬ 
ficients as they act on columns of spectral coefficients from the left. 

Differentiation matrices. 

D'^ : {k,m} ^ {k + l,m + 1} (A.l) 

^t',n = \/2n(n -I- fc -I- TO -I- 1) 5 n’,n-i (A.2) 


D : {fc, to} —>■ {fc -I- 1, m — 1} (A.3) 

= y2(rr+lHMj(r7+r^(5„/_„ (A.4) 

Multiplication matrices. 

R'^ : {k,m} ^ {k,m + 1} (A.5) 


= 


(n -I- TO -I- l)(n -I- fc -I- TO -I- 1) 




•'n'.n ^ (^2n + k + m + l){2n + k + m + 2) ” 

dn'.n—l 


lljl + k) 


{2n + k + m){2n + k + m+ \) 


(A.6) 

(A.7) 


R : {fc, to} —>■ {/c, TO — 1} 


(A.8) 
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i?”/. = 


' {n + l)(n + fc + 1) 

(2n + k -\-m + l)(2n + fc + m + 2) 


^n' ,n+l 


+ A 


{nrn){n-\-k-\-m) 

(2n + k + m){2n + k + m + 1) 


(A.9) 

(A.IO) 


Conversion matrices. 


C : {k, m} —)■ {fc + 1, m} 


Cn'.n — 


' {n + k + l){n +k + m + 1) 

(2n + /c + m + l)(2n + k + m + 2) 

/ n(n + m) 

y (2n + fc + m){2n + k + m+ 1) 


^n' ,n— 1 


(A.ll) 

(A.12) 

(A.13) 


: {k, m} —>■ {fc — 1, m} 


(A.14) 


C^ = 

^n' ,n 


(n + l)(n + TO + 1) 


(2n + k + m+ l){2n + k + m + 2) 


Sn' 


n',n+l 


+ 1 


' (n + A:)(n + A: + m) 

(2n + k + m){2n + k + m + 1) 


(A.15) 

(A.16) 


B. Dirac’s bra-ket notation 

Dirac notation is uncommon in numerical analysis applications outside of quan¬ 
tum mechanics. We believe the notation provides many advantages with regard 
to clarity and abstraction. This short appendix reviews the important aspects 
of the notation as it applies to the bulk of the paper. 

In the appendix, (5„(r) denotes any one of the Jacobi-based orthogonal polyno¬ 
mials in the main section. We drop the k,m indices that parameterise different 
bases. We represent functions of r such that 

OO 

f{r) = (B.l) 

n—0 
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For a given weight function w{r) [e.g.^ r(l — 


Qnir)Qn'{r)w{r)dr = 5n,n'■ 


(B.2) 


Therefore a function’s spectral-expansion coefficients result from projecting against 
a given polynomial 


fn = f /(r)Q„(r)w(r)dr. 
^0 


(B.3) 


We define a (discrete) complete orthonormal basis such that 

OO 

{n\n') = Sn,n', ^ 


(B.4) 


n—0 


The element ( n | denotes an infinite row vector with O’s in the first n entries, a 
1 in the (n -I- l)st entry, and O’s after that. For example. 


(0| = [1,0,0,0,0,0,...], (4| = [0,0,0,0,1,0,...]. 


(B.5) 


In terms of the discrete basis and orthonormal functions, we define the contin¬ 
uous basis element corresponding to evaluation at the point r 


(’”1 = 


n=0 


Therefore, 


Also, 


/ I r )( r |ri;(r)dr = I. 
Jo 


{r\r') = '^Quir)Qn{r') = -^S{r-r'), 


(B.6) 


OO OO 

/ |r)(r|w(r)dr = X! X! II / Qn'{r)Qn'{r)w{r)dr (B.7) 

•^0 n = 0 n '=0 


(B.8) 


(B.9) 


which follows from the Christoffel-Darboux formula for orthogonal polynomials. 
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The orthogonal functions represent the coefficients of the r basis in terms of the 
n basis, 


( r I n) = Qn{r). 


(B.IO) 


Now, 

OO 

f{r) = {r\f) = '^Qn{r){n\f) (B.ll) 

n—0 

where 


fn = {n\f) (B.12) 

represent the coefficients of the state / in the n basis. 


The inner product of functions takes the form 


r*l /•! 

g*(r)f{r)w(r)<ir = / ( 5 | r )( r | / )'u;(r)dr 

3 ^0 


(B.13) 


therefore 


g*{r)f{r)w{r)dr = {g\f). 


(B.14) 
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